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Summary 

This  document  describes  the  rationale  for,  and  preliminary  design  of  a  computer-based 
probabilistic  expert  system  to  predict  the  likelihood  of  burial  and  migration  of  abandoned 
underwater  munitions. 

1.  Introduction 

As  a  legacy  of  years  of  naval  activities,  including  training  and  testing,  there  are  numerous 
current  and  former  Department  of  Defense  (DoD)  aquatic  sites  contaminated  with 
munitions  and  explosives  of  concern  (MEC).  At  many  of  these  inland  water  and  coastal 
areas,  the  risk  of  human  interaction  with  unexploded  ordnance  (UXO)  is  of  grave  concern, 
and  efforts  have  been  mandated  to  manage  and  clean  up  these  contaminated  sites.  The 
Munitions  Response  program  of  the  Strategic  Environmental  Research  and  Development 
Program  (SERDP)  is  tasked  with  developing  innovative  methods  to  remediate  and 
sustainably  administer  areas  polluted  by  discarded  munitions.  Towards  this  end,  SERDP  is 
sponsoring  the  development  of  the  Underwater  Munitions  Expert  System  (UnMES),  for 
predicting  the  location  of  munitions  and  their  degree  of  burial  at  underwater  sites  of 
interest. 

Compared  to  terrestrial  sites,  underwater  environments  are  more  dynamic,  with  UXO1 
often  subject  to  mobility,  burial,  and  excavation  by  bottom  currents  driven  by  waves  and 
tides.  The  extent  of  the  search  area  covered  for  a  region  containing  underwater  UXO  is 
limited  using  present  platform  and  sensor  technology.  Therefore  it  is  advantageous  to  have 
the  ability  to  predict  areas  of  munitions  concentration,  exposure  and  temporal  variability, 
to  support  planning  for  efficient  site  remediation. 

The  predictive  model  UnMES  is  built  on  a  probabilistic  Bayesian  network  whose  construct 
inherently  quantifies  and  tracks  the  uncertainty  of  the  predictions.  This  feeds  naturally  into 
risk  assessment  framework  needed  by  site  managers  for  making  informed  remediation 
decisions.  Toward  this  end,  the  Bayesian  network  in  the  UnMES  is  embedded  in  a 
Geographic  Information  System  (GIS)  framework  to  interpret  and  display  spatial  variation. 


1  The  acronym  UXO  is  used  to  denote  MEC  of  all  types  and  sources. 
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The  choice  of  a  probabilistic  approach  is  appropriate  in  situations  where  there  is  large 
inherent  uncertainty  in  the  predictive  models  and  particularly  in  the  scope  and  range  of 
values  of  the  governing  input  parameters.  The  state  of  the  art  of  deterministic  coastal 
system  models  is  at  present  not  sufficient  to  permit  the  scaling  of  first-principle  physical 
processes,  such  as  turbulence-forced  sediment  suspension  that  occurs  at  the  scale  of  a 
breaking  wave,  up  to  the  larger  coastal  morphological  response  which  occurs  on  timescales 
from  months  to  years  or  even  decades  [Hanson  et  ah,  2003].  In  the  case  of  modeling  UXO- 
contaminated  sites,  there  is  substantial  uncertainty  in  the  boundary  and  forcing  conditions: 
the  initial  population  of  the  munitions  is  usually  not  well-known,  and  the  environment  at 
the  site  over  the  time  since  the  deposition  may  not  have  been  recorded.  Therefore,  we 
must  rely  upon  regional  climatological  distributions  and  assume  most  likely  scenarios. 


The  design  of  UnMES  is  focused  on  providing  quantitative  probabilistic  estimates  that 
address  the  questions  site  managers  need  answered.  For  example,  given  an  estimate  of  the 
initial  distribution  of  UXO,  one  might  want  to  know  from  which  areas  within  the  field  the 
munitions  were  more  likely  to  migrate,  and  in  what  direction,  and  thus  what  specific 
locations  they  might  potentially  aggregate.  Furthermore,  given  knowledge  of  the  types  of 
munitions  onsite,  one  might  want  to  know  which  types  are  more  subject  to  mobilization. 
An  important  prediction  concerns  how  frequently  UXO  are  likely  to  be  buried  versus  proud 
on  the  surface  of  the  seabed,  a  factor  that  significantly  impacts  the  performance  of  sensors 
used  to  detect  and  classify  munitions.  Answers  in  the  form  of  probability  distributions  can 
provide  risk  information,  and  can  be  used  to  evaluate  remediation  alternatives. 


1.1  Scope  of  the  present  study 

The  expert  system  currently  under  development  will  be  a  core  component  in  the  larger 
effort  under  the  SERDP  Munitions  Response  (MR)  Program  to  develop  technologies 
allowing  rapid  assessment  of  large-area  underwater  sites  where  detection  and  remediation 
of  munitions  is  required.  For  successful  implementation,  the  overall  assessment  framework 
will  involve  multiple  components  including: 

1.  Site  specific  environmental  data  and  historical  military  archives. 

2.  UnMES:  probabilistic  estimate  of  UXO  migration  and  burial. 

3.  Detection  sensor  performance  predictions. 

4.  Risk  assessment  based  on  human  exposure  and  removal  issues. 

5.  Decision  support  tool  providing  graphical  mapping  and  visualization  of  UXO 
distributions  and  burial  state  for  site  remediation  managers  and  planners. 
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The  present  effort  is  specifically  on  the  development  of  component  2,  the  UnMES,  whose 
relationship  to  the  larger  system  is  diagramed  in  Figure  1.1.  The  other  components  will  be 
part  of  the  long-term  SERDP/ESTCP  program.  To  demonstrate  and  verify  the  operation  of 
UnMES,  however,  assessments  at  specific  sites  will  be  undertaken,  requiring  representative 
site-specific  data  to  be  assembled  and  utilized. 


Figure  1.1  Conceptual  diagram  of  components  of  the  SERDP  MR  program  required  to 

characterize  risk  and  plan  for  remediation  efforts  at  underwater  munitions  response  sites. 

The  scope  of  the  present  effort  to  develop  of  the  prototype  UnMES  entails  the  following 

tasks: 

1.  Development  of  physical  and/or  data-based  process  models  for  mobility,  burial,  and 
re-exposure  (in  collaboration  with  other  SERDP  MR  studies). 

2.  Development  of  a  probabilistic  BN,  including  rationale  for  the  discretization  of 
parameter  values  within  the  network,  and  methodologies  for  generation  of  the 
conditional  probabilities  that  provide  the  basis  for  UXO  state  predictions. 

3.  Incorporation  of  the  Bayesian  network  into  a  GIS  framework  to  demonstrate  the 
spatial  application  of  probabilistic  prediction. 

4.  Acquisition  of  data  from  representative  sites  of  interest  for  demonstration  and 
verification  purposes  using  extant  databases  and  current  SERDP  studies. 

5.  Demonstration  GIS  mapping  of  UXO  accumulation  and  burial  probabilities  at 
representative  sites. 
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These  development  efforts  are  being  supported  by  modeling,  laboratory  and  field  studies 
under  by  the  SERDP  MR  program,  especially  Friedrichs  (VIMS)  and  Garcia  (U.  Illinois)  for 
Task  1;  along  with  Jenkins  (SIO),  Traykovski  (WHOI),  and  Calantoni  (NRL)  for  Task  4. 


In  order  to  provide  useful  estimates  of  munitions  behavior  probabilities,  sufficient 
environmental  data  is  required  for  understanding  the  geologic  setting  and  primary  wave 
and  current  forcing.  The  more  detailed  and  extensive  the  available  background 
environmental  data  the  more  accurate  and  specific  the  resulting  predicted  probability 
distributions  will  be.  However,  UnMES  is  designed  to  provide  estimated  probabilities  even 
when  the  input  description  is  broad  or  generic. 


2  Design  of  UnMES 


The  central  paradigm  of  UnMES  is  a  Bayesian  Network  [Pearl,  1988],  a  graphical 
probability  model  that  reflects  current  knowledge  of  the  likely  behavior  of  UXO  migration 
and  burial  in  a  given  environment,  as  illustrated  in  Figure  2.1.  This  core  construct  should 
be  placed  in  a  2-D  map  framework  in  order  to  accommodate  spatially  varying  input 
parameters,  and  to  distribute  resulting  predictions  in  a  map-based  management  tool. 
Populating  the  2-D  framework  with  the  relevant  environmental  data  requires  building 
connections  from  UnMES  to  a  Geographic  Information  System  (GIS)  developed  for  the  site 
of  interest.  An  overview  of  this  software  framework  is  diagrammed  in  Figure  2.2.  The 
prototype  system  now  under  development  will  be  exercised  at  sites  where  input  data  is 
currently  available.  Development  of  a  full-featured  management  tool  where  in  the  UnMES 
in  embedded  as  a  module  with  GIS  support  [e.g.  Holland,  2014]  remains  a  task  for  a  future 
phase  of  the  SERDP  MR  program. 

The  focus  of  the  preliminary  UnMES  design  will  be  on  wave-dominated  climates,  where  the 
important  environmental  forcing  is  primarily  due  to  the  bottom  orbital  flows  from  wind- 
driven  waves  -  a  situation  of  primary  interest  for  remediation  efforts;  and  an  environment 
in  which  several  SERDP  field  experiments  have  been  conducted  and  can  provide  validation 
data  [Calantoni  et  ah,  2015,  Traykovski,  2015].  In  these  regions,  encompassing  the  coastal 
ocean  and  large  lakes,  the  sediment  is  generally  comprised  of  non-cohesive  sands. 
Subsequent  to  the  development  and  proofing  of  this  prototype  version,  the  design 
framework  can  be  expanded  to  include  the  other  underwater  environments  of  interest, 
including  estuaries  and  rivers,  which  have  different  circulation  patterns  and  include 
additional  forcing  mechanisms,  such  as  tides,  flood  erosion  and  fluvial  deposition.  Sand 
sediments  under  waves  are  generally  resistant  to  penetration,  whereas  in  estuaries  and 
rivers  the  bottom  type  is  often  soft  clays  or  muds  where  burial  upon  initial  impact  with  the 
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seabed  may  be  important.  Additional  wave-driven  regions  with  different  underlying 
geology,  such  as  large  cobbles,  and  coral  reefs,  can  be  added  to  future  versions  of  UnMES. 


For  this  preliminary  UnMES  design,  the  focus  is  on  seasonal  and  annual  timescales  where 
the  process  probabilities  can  be  treated  statistically.  Modeling  the  dynamic  evolution  in 
time  steps  over  a  specified  time  horizon  will  be  explored  in  a  more  advanced  version. 


Figure  2.1.  Diagram  of  the  core  Bayesian  Network  (BN)  in  UnMES  modeling  the  burial  and 
migration  response  in  a  given  environmental  setting  (specified  by  blue  nodes).  Knowledge  of 
the  quantity,  type,  and  initial  burial  condition  of  the  munitions  are  provided  in  the  orange 
nodes.  The  resulting  predictions  are  found  in  the  output  (green)  nodes. 

2.1  Introduction  to  Bayesian  Networks 

The  core  construct  of  UnMES,  the  Bayesian  Network  (BN),  sometimes  called  a  Bayesian 
belief  network  [Pearl,  1988],  is  formed  as  a  directed  acyclic  graph  (DAG)  of  variables 
(nodes)  connected  by  directed  links  (arrows).  Each  node  represents  a  random  variable  and 
the  links  represent  statistical  dependencies  between  them.  As  the  value  at  each  node  is  a 
probability  mass  function,  a  BN  is  a  method  to  model  a  multivariate  domain  that  contains 
uncertainty.  This  uncertainty  can  be  due  to  incomplete  specification  of  the  interactions 
governing  the  behavior  of  the  domain,  such  as  the  simplification  of  wave  orbital  velocities 
using  linear  theory,  or  incomplete  knowledge  of  the  state  of  some  of  the  variables,  as  when 
water  depth  is  estimated  from  an  old  or  poorly-sampled  bathymetry.  The  uncertainty  is 
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propagated  through  the  BN  by  the  probabilistic  inference  using  Bayes  theorem2.  Figure  2.1 
shows  the  DAG  of  nodes  that  form  the  core  BN  of  UnMES  for  a  wave-dominated 
environment. 


The  BN  portion  of  UnMES  is  designed  and  implemented  using  the  commercial  software 
package  'Netica'  [Norsys,  2015],  which  provides  a  convenient  user  interface  for  graphically 
designing,  building  and  modifying  networks,  as  well  as  entering  new  input  probabilities 
("findings”)  and  displaying  the  resulting  predictions.  To  represent  the  spatial  variability 
across  a  contaminated  underwater  site,  the  core  Bayesian  Network  is  replicated  at  cells 
across  a  two-dimensional  representation  of  the  geographic  region  of  interest  as 
diagrammed  in  Figure  2.2.  The  overall  structure  is  implemented  in  Matlab  [MathWorks 
Inc.,  2014],  which  interacts  with  Netica  through  its  JAVA  Application  Programmer  Interface 
(API).  Further  details  are  provided  in  Sections  2.2  and  5. 


Bathymetry 


GIS  -  Input  Properties 


Wave  Height 

UXO  tvPe  and  abundan 
etc ... 

Province  into  cells 
with  unique  input 
combinations 


Replicate  UnMES 
BN  at  each  cell 


Figure  2.2  Conceptual  illustration  of  GIS  software  framework  for  spatial  implementation  of 
UnMES. 


After  a  Bayesian  Network  has  been  constructed,  it  can  be  applied  to  a  particular  situation 
by  entering  information  (known  as  "evidence"  or  "findings"),  setting  the  values  of  nodes 
that  are  known.  For  example,  in  UnMES,  these  nodes  will  include  the  environmental 
observations  for  a  particular  location.  Then  probabilistic  inference  (also  called  "belief 
updating"  in  a  BN)  determines  new  probabilities  for  the  states  of  all  the  other  dependent 


1  For  background  on  Bayesian  reasoning  and  its  application  to  graphical  models  see  Jensen  and  Nielsen  [2007]. 
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nodes.  These  are  displayed  in  Netica  as  "belief  bars"  which  show  the  shape  of  the 
probability  distribution.  Before  performing  belief  updating,  the  BN  must  be  compiled.  For 
this,  Netica  uses  an  efficient  junction  tree  algorithm  based  the  assumption  of  conditional 
independence  [Roller  and  Friedman,  2009].  In  order  to  implement  the  inference  algorithm 
in  Netica,  every  node  must  be  discrete  with  a  finite  number  of  states.  Some  nodes  are 
naturally  discrete,  e.g.  UXO  type.  However  most  environmental  variables  are  inherently 
continuous,  so  that  total  range  of  the  variable  must  be  broken  into  a  number  of  intervals,  or 
states,  whose  sizes  do  not  need  to  be  constant.  The  intervals  should  be  small  enough  to 
discriminate  different  response  states,  but  not  so  narrow  that  the  intervals  have  a  higher 
resolution  than  the  available  data  warrant  and/or  such  thatthe  computational  requirement 
is  prohibitive.  The  discretization  choice  for  each  node  is  an  important  part  of  the  design 
process  and  will  be  discussed  in  more  detail  in  Section  3  for  the  present  UnMES  application. 


2.2  Nodes,  Links  and  Conditional  Probability  Tables 


A  node  in  a  BN  with  no  links  coming  into  it  is  considered  a  "parent"  node,  and  will  contain 
an  unconditional  probability  table  quantifying  the  probability  mass  function  (PMF)  over 
the  states  of  the  variable  that  it  represents.  In  the  UnMES,  these  parent  nodes  act  as  input 
variables,  providing  the  boundary  conditions  for  an  environmental  scenario,  for  example, 
the  Wave  Height  and  UXO  Type  nodes  shown  in  Figure  2.1.  If  a  node  has  one  or  more 
links  pointing  into  it,  it  is  a  "child”  node,  and  has  a  conditional  probability  table  (CPT).  Each 
entry  in  the  CPT  contains  a  conditional  probability  for  a  specific  state  of  that  child  node, 
given  a  specific  configuration  of  the  states  of  its  parents.  In  UnMES,  there  are  also 
intermediate  variables  such  as  the  computed  Bottom  Orbital  Velocity,  which  is  a  child  to 
the  wave  forcing  nodes,  and  also  a  parent  to  the  node  ScourBurial. 


A  Bayesian  Network  can  be  used  in  diagnostic  mode,  i.e.  findings  can  be  specified  at  a 
"results”  node  (e.g.  migration  distance),  and  then  the  most  probable  configuration  of  values 
for  the  rest  of  the  network  can  be  determined.  Another  important  diagnostic  feature  that 
Netica  provides  is  the  ability  to  do  sensitivity  analysis,  producing  a  report  of  how  much  the 
distribution  of  states  at  a  selected  node  could  be  influenced  by  a  particular  specified  state 
at  each  of  the  other  nodes  in  the  net.  The  use  of  sensitivity  analysis  will  be  presented  in  the 
forthcoming  Demonstration  Report  on  the  prototype  UnMES  implementation  [Rennie  and 
Brandt,  2015b].  Both  this  single-finding  sensitivity  analysis,  and  a  comparison  of  results 
from  inputs  of  varying  levels  of  detail  can  identify  which  factors  contribute  the  most  to  the 
predicted  uncertainty.  Because  environmental  data  acquisition  is  often  the  most  expensive 
part  of  a  site  remediation  study,  this  evaluation,  performed  early  in  the  process,  would  be 
important  to  help  refine  data  requirements. 
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Relationships  in  the  network  are  defined  by  each  node’s  CPT,  which  is  a  multi-dimensional 
array  that  has  one  probability  value  for  each  combination  of  parent  and  child  states.  The 
CPT  will  be  an  N-dimensional  table,  where  N  is  the  number  of  parents,  which  has  a  a  total 
number  of  entries  equal  to  MrlVh-  ...  Mn,  where  Mi  is  the  number  of  states  in  parent  node  i. 
The  number  of  entries  in  the  CPT  will  be  the  multiple  of  the  number  of  child  states  times 
the  number  of  states  in  each  parent.  Clearly,  if  variables  are  discretized  into  numerous 
intervals  (states),  the  size  of  the  CPT  grows  rapidly. 

Most  of  the  variables  in  UnMES  describe  physical  processes  and  are  continuous  in  nature. 
The  choice  of  the  discrete  intervals  by  which  to  represent  each  node  should  correctly 
reflect  the  sensitivity  of  its  response  to  its  forcing,  yet  cover  the  range  of  interest  in  as  few 
states  as  possible,  so  as  to  be  represented  by  a  compact  CPT  [Plant  and  Holland,  2011]. 
Note  that  there  is  no  requirement  that  the  intervals  be  the  same  size  across  the  range  of  the 
variable.  When  designing  UnMES,  the  resolution  of  the  discretization  interval  reflects  a 
trade-off  between  resolving  important  information  and  avoiding  over-specification  of 
details  occurring  at  scales  smaller  than  the  expected  uncertainty.  The  choice  of  state 
intervals  is  an  important  design  decision;  the  rationale  for  UnMES  state  boundaries  is 
detailed  in  Section  3. 

Populating  each  of  the  entries  in  the  CPT  is  the  main  effort  in  building  the  BN.  Often  CPT 
population  in  Bayesian  Networks  is  accomplished  by  training  on  multiple  data 
observations.  An  example  application  of  BN  to  predict  coastal  response  to  sea  level  rise  by 
Gutierrez  et  al.  [2011]  was  trained  on  a  large  historical  data  set  developed  by  the  U.S. 
Geological  Survey  (USGS).  As  there  are  only  a  modest  amount  of  observations  relating  to 
the  migration  and  burial  of  underwater  munitions,  the  CPTs  in  UnMES  will  be  largely  filled 
in  by  a  Monte  Carlo  exploration  over  the  domain  of  interest  using  simple  deterministic 
process  models.  The  deterministic  models  are  "simple",  in  that  they  are  focused  on 
capturing  the  first-order  response,  and  do  not  attempt  to  replicate  small  details  of  the 
behavior,  given  our  insufficient  knowledge  of  the  system.  Available  laboratory  and  field 
data  is  used  to  calibrate  the  parameters  of  these  models  [Rennie  and  Brandt,  2014].  The 
state  of  the  deterministic  models  developed  for  use  in  the  UnMES  is  described  in  Section  4. 

2.3  Spatial  Distribution:  Replicas  of  the  Bayesian  Network 

At  UXO  contaminated  sites,  managers  may  need  to  answer  the  question  of  whether  certain 
areas  within  the  site  are  more  susceptible  to  munitions  migration/aggregation  or  burial 
than  others;  or  more  specifically,  what  are  the  probabilities  of  UXO  aggregation  and  burial 
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for  local  areas  within  the  site  of  interest.  Another  prediction  requiring  spatial  mapping 
occurs  when  regions  of  higher  munitions  contamination  (i.e.  larger  abundance  of  UXO)  are 
identified  in  the  initial  site  conceptual  model,  or  if  there  are  sub-regions  where  the 
potential  for  human  interaction  is  particularly  great.  To  address  these  questions  it  is 
important  to  implement  the  UnMES  in  a  2-D  spatial  configuration  that  will  naturally 
connect  to  map-oriented  risk  assessment  products. 


The  application  of  Bayesian  networks  in  a  GIS  framework  for  geospatial  reasoning  has  been 
explored  in  several  fields,  including  hydrology  [USEPA,  2007],  ecosystem  modeling 
[Landuyt,  et  al,  2015],  and  military  logistics  [Laskey,  et  al.,  2010].  However,  the  experience 
with  implementation  of  this  construct  is  relatively  new,  with  no  established  software 
framework.  Norsys  is  in  the  process  of  developing  a  product,  GeoNetica,  which  facilitates 
the  use  of  Netica  Bayesian  Networks  in  conjunction  with  raster-style  GIS  input  and  output. 
GeoNetica  interfaces  with  multiple  commercially  available  GIS  systems  and  will  be  of  great 
benefit  for  implementing  the  UnMES  in  coordination  with  the  site-specific  environmental 
GIS  databases  (Figure  2.2).  GeoNetica  is  in  "Beta"  release  to  a  limited  customer  base. 
Recently,  JHU/APL  was  approved  to  join  this  initial  test  group.  For  this  preliminary 
implement  of  UnMES,  we  will  couple  the  UnMES  Netica  BN  with  the  open-source 
Geographic  Information  System  QGIS  [www.qgis.org]  to  demonstrate  the  spatial 
framework  analysis  approach.  Practical  application  of  UnMES  would  require  its  integration 
into  a  GIS  framework  that  meets  the  geospatial  data  standards  required  by  the  agency  with 
cognizant  authority  for  the  contaminated  property. 


As  with  the  continuous  nodes  in  the  BN,  the  spatial  domain  should  be  discretized  into 
regions  of  relative  uniformity,  analogous  with  the  nodes’  state  intervals.  As  with  the  state 
intervals,  the  uniform  geographic  regions  are  not  necessarily  of  all  of  the  same  size.  Often 
spatial  areas  are  partitioned  into  a  uniform  grid  of  evenly  spaced  regions  (pixels),  as  in 
raster  GIS  format.  This  is  a  straightforward  approach  which  is  simple  to  work  with,  but  can 
be  inefficient,  since  large  numbers  of  pixels  may  contain  the  (approximately)  same  feature 
value.  The  environmental  data  required  by  UnMES,  such  as  bathymetry,  is  frequently 
available  as  a  GIS  raster  layer,  and  should  be  preprocessed  to  a  pixel  resolution  no  finer 
than  needed  to  resolve  the  state  intervals  in  the  corresponding  BN  node.  Wave  information 
from  gridded  coastal  models  is  naturally  stored  in  raster  (pixel)  form. 


Preprocessing  steps  are  necessary  to  prepare  the  GIS  data  before  the  input  nodes  of  the  BN 
can  be  populated.  The  several  input  layers  must  be  reshaped  for  common  pixel  alignment 
and  resolution.  There  are  two  approaches  to  implementing  a  BN  in  a  GIS  framework:  a  fast 
approach  which  is  memory-intensive,  or  a  slow  one  with  less  memory  usage  [Landuyt,  et 
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al.,  2015].  The  slow  mode  is  based  on  GeoNetica  and  runs  the  BN  model  at  each  pixel 
independently.  For  the  fast  approach,  collections  of  pixels  which  represent  unique 
combinations  will  be  combined  and  referred  to  as  a  "cell".  There  usually  will  be  many  fewer 
cells  than  pixels  for  a  given  site  implementation.  Then  the  BN  is  replicated  at  each  of  these 
smaller  set  of  locations.  This  approach  requires  a  catalog  that  lists  all  unique  cells  with 
their  corresponding  pixels,  and  then  maps  the  probabilistic  model  output  back  into  raster 
form.  This  catalog  allocates  considerable  amounts  of  memory.  However,  the  unique  cell 
approach  allows  for  customizing  the  nodes’  discretization  intervals  in  the  nodes  for  more 
efficient  dynamic  range  coverage  (Section  3).  The  advantages  of  each  approach  will  be 
evaluated  as  we  gain  experience  with  implementation  of  UnMES  at  additional  sites. 


In  the  initial  version  of  UnMES,  the  core  BN  is  run  for  each  cell  of  the  joint  GIS  database 
independently.  The  probabilistic  output  is  redistributed  back  through  the  catalog  to  the 
corresponding  pixel  locations  in  output  GIS  raster  layers  that  can  be  displayed  as  maps 
visualizing  the  predictions  and  their  corresponding  uncertainty.  The  presentation  of  results 
can  be  displayed  in  several  forms,  including  maps  of  most  probable  state  (or  mode  state), 
accompanied  by  a  map  of  the  probability  of  that  reported  state;  or  a  map  of  expected  value 
along  with  a  second  map  of  standard  deviation.  A  discussion  of  the  methods  of  visualizing 
uncertainty  is  presented  in  Section  5. 


2.4  Temporal  Evolution 


The  process  of  UXO  migration  on  the  seabed  is  envisioned  to  be  a  repeated  sequence  of 
movement,  then  burial,  followed  eventually  by  re-exposure,  leading  to  possible  further 
movement.  The  ESTCP  field  tests  conducted  at  the  Field  Research  Facility  (FRF)  in  Duck, 
NC  [Wilson  et  al,  2008a]  observed  surrogate  UXO  buried  under  sediment  to  depths  several 
times  the  object  diameter.  Subsequent  observations  revealed  these  UXO  were  at  different 
locations,  but  again  deeply  buried.  It  is  inferred  that  local  bathymetric  change,  either  from 
bedform  migration  or  shoreface  adjustment,  caused  the  UXO  to  become  unburied,  allowing 
the  UXO  to  migrate,  and  then  rebury.  Based  on  a  compilation  of  field  test  results  from 
SERDP  researchers  [Calantoni  et  al,  2015,  Traykovski,  2015]  as  well  as  modeling  studies 
by  jenkins  et  al  [2013],  burial  is  by  far  the  most  probable  state  of  the  UXO  at  any  given 
time,  with  re-exposure  and  movement  occurring  very  infrequently  at  most  sites.  However, 
it  may  be  just  these  unusual  events  that  are  of  most  concern  to  site  managers.  At  many  of 
the  sites,  the  munitions  have  been  in  residence  for  years,  if  not  decades.  The  potential  for 
infrequent  events  to  occur,  though  small,  is  finite. 


10 


Johns  Hopkins 

APPLIED  PHYSICS  LABORATORY 


FPS-T- 15-0333 
Technical  Report  Page  1 1 


The  approach  taken  by  UnMES  is  to  model  the  burial-migration  cycle,  not  explicitly  as  a 
time-stepped  sequence  of  processes,  but  as  the  combined  probability  of  burial  and  motion 
states,  given  a  long-term  distribution  of  environmental  forcings.  We  focus  here  on  a  wave- 
forced  environment,  so  the  temporal  distributions  of  wave  heights  and  their  corresponding 
periods  are  required  input,  along  with  a  quantification  of  the  probability  of  changes  in 
seabed  elevation  that  provide  the  potential  for  re-exposure.  The  time  scale  of  temporal 
evolution  is  implicit  in  the  duration  over  which  the  statistics  are  derived.  For  example,  if 
the  distribution  of  seabed  elevation  changes  are  computed  from  measurements  or  modeled 
processes  representing  weeks  or  a  few  months,  that  PDF  captures  re-exposure  potential 
due  to  local  bedform  migration,  but  not  that  caused  seasonal  shoreface  adjustment,  which 
occurs  over  annual  time  scales. 


3  Details  of  Nodes  in  UnMES 

The  nodes  in  the  core  UnMES  BN  can  be  grouped  into  three  sets.  In  a  predictive  mode, 
normally  the  input  nodes  which  specify  the  geologic  and  meteorological  setting  (shown  as 
blue  in  Figure  2.1),  will  have  their  states  set  using  values  from  the  GIS  environmental 
database.  There  are  several  variables,  denoted  by  orange,  that  provide  existing  knowledge 
of  the  initial  disposition  of  the  UXO  that  are  contaminating  the  site,  such  as  the  kinds  of 
munitions  and  their  relative  abundance. 


The  probability  distributions  resulting  in  the  output  prediction  nodes  (green  set  shown  in 
Figure  2.1)  contain  the  information  that  will  feed  into  site  manager  risk  assessment  tools. 
The  set  of  intermediate  nodes  is  where  the  functionality  of  UnMES  lies.  The  CPTs  for  these 
nodes  are  filled  in  using  our  knowledge  of  ocean  physics  along  with  recent  laboratory  and 
field  research  into  the  behavior  of  munitions  on  the  seabed  in  the  form  of  process  models 
(see  Section  4)  and  rule  based  expert  input  for  situations  where  physics  based  models  are 
not  available. 


A  crucial  step  in  constructing  the  BN  is  the  choice  of  discretization  intervals  for  both  the 
input  and  intermediate  nodes.  The  primary  determinant  of  the  discretization  choice  is  the 
need  to  adequately  resolve  the  behavior  of  the  UXO.  However,  in  order  to  maintain  a 
reasonable  size  for  the  multi-dimensional  CPTs,  there  is  a  need  to  limit  the  number  of 
discretization  intervals.  For  sufficient  discrimination  of  the  munitions  burial  and  migration 
patterns,  a  large  number  of  states  may  be  required,  but  the  inherent  uncertainty  of  input 
measurements  should  be  considered  before  creating  unrealistically  high-resolution 
intervals.  Whitehouse  [1998]  cites  typical  uncertainties,  ranging  from  10%  to  20%  for  the 
input  environmental  parameters  in  UnMES,  when  there  has  been  an  in-situ  field  survey.  So 
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the  uncertainty  when  using  generic  regional  predictions  can  be  expected  to  be  significantly 
larger. 

The  precision  with  which  UnMES  needs  to  predict  burial  is  determined  by  a  number  of 
factors.  Doctrinal  burial  categories  used  in  the  Navy’s  Mine  Warfare  and  Environmental 
Decision  Aids  Library  (MEDAL),  which  reflect  the  sensitivity  of  various  underwater  mine¬ 
hunting  sensors,  are  set  to  intervals  with  boundaries  of  10%,  20%,  75%  and  100%  burial. 
In  addition  to  supporting  sensor  performance  prediction,  UnMES  also  calculates  when  UXO 
are  buried  sufficiently  so  that  mobility  is  suppressed.  The  percentage  burial  at  which  this 
"lockdown"  occurs  has  been  estimated  to  be  between  40%  to  50%  [Wilson  et  al,  2008a]. 
Therefore  the  scour  burial  states  of  10%,  20%,  50%,  80%  and  100%  of  the  UXO  diameter 
are  proposed  for  the  initial  version  of  UnMES.  These  discretization  intervals  are  denoted  on 
Figure  3.1,  and  listed  in  column  1  of  Table  3.1. 


Table  3.1  Discretization  Intervals  for  some  UnMES  Bayesian  Network  Nodes 


Scour 

Burial 

(%) 

Bottom 

Velocity  Ubot 
(cm/s) 

Water 
Depth  h 
(m) 

Wave 
Height  (m) 
for  /7=2.5m 

Wave 
Period  (s) 
for  b=2.5m 

Wave 
Height  (m) 
for  h  =  7  m 

Wave 
Period  (s) 
for  h  =  7  m 

<  10 

<  10 

1  to  1.5m 

<0.3 

<4 

<0.4 

<4 

10  to  20 

10  to  20 

1.5  to  2m 

0.3  to  0.4 

>4 

0.4  to  0.7 

4  to  6 

20  to  50 

20  to  30 

2  to  2.5m 

0.4  to  0.5 

0.7  to  1.0 

>6 

50  to  80 

30  to  40 

2.5  to  3m 

0.5  to  0.65 

1.0  to  1.4 

80  to  100 

40  to  50 

3  to  4m 

0.65  to  0.8 

1.4  to  1.7 

fully 

50  to  60 

4  to  5m 

0.8  to  1.0 

1.7  to  2.0 

buried 

60  to  75 

5  to  6m 

1.0  to  1.2 

2.0  to  2.5 

75  to  90 

6  to  8m 

1.2  to  1.4 

2.5  to  3.0 

90  to  105 

8  to  10m 

1.4  to  1.7 

3.0  to  3.5 

105  to  120 

10  to  15m 

1.7  to  2.0 

3.5  to  4.0 

120  to  160 

2.0  to  2.6 

4.0  to  5.3 

surf  zone 

>  2.6m 

>  5.3  m 

The  functional  response  to  forcing  is  highly  sensitive  for  many  of  the  relationships 
modeled  here,  but  the  domain  of  strong  sensitivity  to  one  variable  can  change  depending 
on  the  value  of  a  second  variable.  This  will  be  illustrated  in  detail  for  the  Bottom  Velocity 
node  in  the  following  section.  Because  the  UnMES  GIS  framework  creates  a  new  instance  of 
the  BN  at  each  spatial  cell,  we  can  take  advantage  of  the  reduction  in  dynamic  range  for 
variables  that  are  constrained  within  that  spatial  region,  such  as  water  depth  and  sediment 
grain  size.  For  example,  because  waves  break  in  water  depths  on  the  order  of  their  height, 
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the  Wave  Height  node  for  the  BN  replica  at  spatial  cells  in  very  shallow  water  need  not 
include  states  for  the  larger  wave  heights  that  may  be  found  farther  offshore. 


In  a  previous  implementation  of  a  Bayesian  network  to  predict  mine  burial  [Rennie  et  al, 
2007]  the  CPT  was  defined  using  an  off-line  Monte-Carlo  computation  of  a  predictive 
deterministic  model.  For  UnMES,  several  of  the  deterministic  models  (see  Section  4)  are 
simple  enough  to  specify  the  equation  defining  the  relationship  directly  in  Netica’s  node 
Equation  box  whereby  Netica’s  Equation-to-Table  option  can  be  used.  The  Equation-to- 
Table  option  performs  a  similar  Monte-Carlo  sampling  of  all  states  on-the-fly  to  build  the 
table  internally.  This  makes  it  possible  to  construct  each  BN  replica  in  a  custom  manner  for 
each  spatial  cell  in  the  joint  GIS  database  (See  Section  3.3.2). 

3.1  Rationale  for  Discretization:  Bottom  Velocity 

The  bottom  velocity  (Ubot)  is  an  intermediate  node  that  acts  as  a  parent  to  both  the  burial 
and  migration  nodes.  Velocity  at  the  bottom  of  the  water  column  is  the  peak  orbital  velocity 
computed  from  linear  wave  equations  [Soulsby,  1997],  as  a  function  of  the  value  of  the 
input  states:  Wave  Height,  Period  and  Water  Depth.  To  determine  the  appropriate  choice  of 
state  interval  width  for  the  Bottom  Velocity  node,  the  sensitivity  of  the  parameterized 
model  for  scour  burial  is  considered.  This  model  was  the  focus  of  recent  studies  [Rennie 
and  Brandt,  2014],  reviewed  in  Section  4.1. 

Across  a  given  spatial  cell,  the  sediment  characteristics  are  fixed  to  a  single  value.  In  Figure 
3.1,  scour  burial,  %B  (in  percent)  is  plotted  against  increasing  bottom  velocity  for  two 
sediment  values,  and  for  two  different  UXO  shapes.  The  type  "bomb"  is  used  to  represent  a 
larger  cylindrical  UXO,  nominal  diameter  of  150  mm,  and  "shell"  denotes  artillery  with  a 
tapered  shape  of  nominal  diameter  80  mm.  The  model  for  %B  is  dependent  on  the 
parameter  0,  the  sediment  Shields  number  [Friedrichs,  2014],  which  is  a  function  of  the 
sediment  grain  size  and  Ubot2.  The  "bomb"  curves  show  the  relationship  %B  =  100-2. 2-00-85, 
the  best  fit  parameters  found  from  the  data  on  cylinder  burial  in  Rennie  and  Brandt  [2014]. 
The  two  sediment  sizes  bound  most  of  the  range  of  sand  grains,  from  dsed  =  0.2  mm  (very 
fine  sand)  to  dsed  =  0.8  mm  (coarse  sand).  The  "shell”  curves  show  %B  =  100-5-0,  the 
empirical  fit  based  on  the  burial  of  tapered  shapes.  At  cells  with  the  very  small  grain  size, 
burial  is  highly  sensitive  to  bottom  velocity,  especially  for  the  "shell"  burial  with  the  larger 
scale  factor  applied.  However,  for  a  "bomb"  in  coarse  sand,  it  takes  a  sizable  increase  in 
bottom  velocity  to  cause  substantial  burial. 
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To  understand  burial  behavior  to  the  precision  proposed  in  Table  3.1,  Ubot  should  be 
resolved  to  within  approximately  0.05  m/s  intervals  for  a  cell  with  a  sediment  type  of  very 
fine  sand  (red  and  black  lines  in  Figure  3.1).  For  a  BN  with  dsed  set  to  coarse  sand,  the 
velocity  intervals  could  be  as  wide  as  0.15  m/s.  In  addition  to  considering  the  resolution 
required  to  discern  changes  within  the  deterministic  relation,  the  uncertainty  of  the 
empirical  fit  (Section  4.1)  should  be  considered.  Overlaid  on  Figure  3.1  are  the  95% 
confidence  ranges  (dashed  lines)  for  the  empirical  fits  for  burial  behavior  of  the  two  UXO 
shapes.  Because  there  were  very  few  sample  data  points  for  the  tapered  shapes,  the 
uncertainty  is  noticeable  larger  for  the  "shell"  relationship.  Taking  the  95%  confidence 
ranges  into  account,  a  velocity  interval  of  0.10  m/s  accurately  reflects  current  knowledge  of 
the  percent  burial  relationship  with  Ubot.  These  Bottom  Velocity  intervals  are  indicated  by 
marks  along  the  left  half  of  the  x  axis  in  Figure  3.1. 


Figure  3.1  Scour  burial  as  a  function  of  bottom  velocity  for  a  bomb-like  object,  and  for  a 
tapered  shell.  The  green  and  blue  lines  show  burial  in  coarse  sand,  and  the  red  and  black 
lines  show  burial  in  very  fine  sand.  The  dashed  lines  indicate  95%  confidence  about  the 
slope  parameter  from  empirical  fits.  The  grey  horizontal  bars  delimit  the  scour  burial  state 
boundaries  at  10%,  20%,  50%,  80%  and  100%.  Grey  triangles  mark  Ubot  state  intervals. 

For  velocities  >1.2  m/s,  the  slowly-burying  munitions  will  become  fully  buried  even  in 
coarse  sands.  This  discrimination  interval  of  larger  Ubot  is  of  interest  for  estimating 
migration  (Section  3.7)  in  the  situation  where  erosional  processes  have  unearthed  the  UXO. 
Because  our  knowledge  of  migration  distances  is  limited,  the  choice  of  wide  Ubot  intervals  at 
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the  higher  speeds  is  appropriate.  With  bottom  orbital  velocities  above  about  2  m/s,  the 
wave  starts  to  break,  and  linear  wave  formulation  is  no  longer  applicable.  The  physics  of 
water  motion  in  the  surf  zone  is  extremely  complex,  and  for  the  present  we  will  simply 
assign  a  "high"  value  to  Ubot  for  these  conditions.  For  the  initial  version  of  UnMES,  we  set 
the  state  interval  widths  for  Ubot  at  0.1  m/s  up  to  0.6  m/s,  increasing  to  0.15  m/s  widths  up 
to  1.2  m/s,  followed  by  0.4  m/s  interval  to  1.6  m/s.  All  larger  values  will  be  denoted  "surf 
zone"  and  assigned  Ubot  >  2.0  m/s,  as  shown  in  Table  3.1.  These  state  intervals  are 
illustrated  in  Figure  3.2  where  the  nodes  are  expanded  to  show  their  "belief  bars",  the 
method  that  Netica  uses  to  display  the  updated  PMF  at  each  node.  Note  that  the  Wave 
Height  distribution  represents  the  temporal  variation  within  a  given  GeoNetica  cell,  while 
Water  Depth  remains  fairly  constant  within  the  cell.  Further  discussion  of  Figure  3.2  is 
presented  in  Section  5. 
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Figure  3.2  Netica  BN  for  scour  burial.  Representative  beliefs  bars  of  the  BN  nodes  for  waves, 
water  depth,  bottom  orbital  velocity  and  percent  burial  by  scour,  showing  the  discretization 
choices  defining  the  state  intervals.  Water  depth  is  set  to  100%  belief  for  the  BN  replica 
within  a  spatial  cell  defined  by  isobathic  depths  between  2.5  and  3  m.  The  UXO  Type  node  is 
set  to  represent  an  area  where  all  underwater  munitions  are  similar  to  naval  artillery  shells. 
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3.2  Input  Nodes:  Geologic  Setting 

The  geological  characteristics  of  the  location  are  specified  in  the  UnMES  by  a  node 
specifying  the  sediment  grain  size,  dsed,  along  with  a  bathymetric  variable,  the  node 
WaterDepth.  A  third  node,  called  hSigma,  quantifies  the  vertical  variability  in  bathymetry 
about  the  mean  water  depth,  h.  These  variables  are  assumed  to  vary  over  the  spatial 
domain,  but  are  temporally  static.  Sediment  grain  size  will  be  constant  within  a  given  cell, 
while  the  temporal  variability  of  h  is  captured  in  the  hsigma  statistic.  The  map  region  to  be 
modeled  is  partitioned  into  cells  based  on  the  merged  GIS  raster  layers  for  these  variables, 
as  discussed  in  Section  2.2  so  that  in  each  cell,  these  nodes  have  a  fixed  value;  i.e.  the  input 
PMF  is  100%  in  a  single  state.  By  narrowing  the  range  of  several  input  (parent)  nodes,  the 
size  of  the  CPTs  of  the  intermediate  (child)  nodes  is  reduced. 

3.2.1  Sediment  characterization 

The  characteristic  of  the  sediment  is  given  by  a  single  parameter,  dsed,  the  diameter  of  the 
sand  grains  in  mm.  At  sites  of  interest,  a  sediment  analysis  would  usually  be  performed, 
where  separation  by  sieving  results  in  the  PMF  of  sand  grain  size  (e.g.  Wilson  et  al.,  2008a). 
Because  the  scour  process  is  sensitive  to  sediment  size  (as  seen  in  Figure  3.1),  the  details  of 
the  grain  size  distribution  were  shown  to  play  a  role  in  CFD  studies  of  UXO  burial  [Jenkins 
et  al.,  2012].  The  UnMES  approach  is  focused  on  first-order  effects,  under  the  assumption 
that  the  overall  uncertainty  (due  to  imperfect  knowledge  of  the  location  and  history) 
renders  details  irrelevant.  Therefore  the  grain  size  distribution  is  represented  by  its 
median  value.  This  value  is  treated  by  Netica  as  a  constant,  rather  than  a  state  with  a  given 
interval  width,  so  no  discretization  choice  is  required  for  dsed ■  Note  that  nodes  in  Netica 
which  represent  constants  do  not  require  links,  but  are  available  as  input  to  all  other  nodes. 

For  many  U.S.  coastal  areas,  there  is  sediment  database  available  from  NAVOCEANO.  A  map 
of  the  NAVOCEANO  Bottom  Sediment  Type  data,  forming  a  QGIS  raster  layer,  is  shown  in 
Figure  3.3  for  the  vicinity  of  Panama  City,  FL.  In  many  wave-dominated  locations,  the  sand 
will  be  fairly  well  sorted  by  the  wave  action,  so  that  water  depth  and  sand  grain  size  are 
inversely  correlated.  However,  due  to  geologic  history,  outcroppings  of  different  sediments 
can  occur,  as  seen  in  Figure  3.3,  where  pockets  of  gravel  and  shell  are  interspersed. 
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Figure  3.3  QGIS  map  of  bottom  sediment  type  from  NAVOCEANO  database  from  the  region  in 
the  Gulf  of  Mexico  near  Panama  City,  Florida.  The  sediment  type  can  be  mapped  to  a  median 
sand  grain  size,  which  is  assumed  to  not  vary  with  time.  Red  circles  mark  the  locations  of  the 
TREX13  field  experiment  locations  [Calantoni,  2015]. 


3.2.2  Water  Depth 

Coastal  bathymetry  is  generally  available  at  all  U.S  sites  at  a  minimum  as  charts  from 
NOAA’s  Office  of  Coastal  Survey  and  digitally  from  the  from  the  National  Ocean  Service 
(NOS)  hydrographic  survey  archive.  Another  source  is  the  U.S.  Naval  Oceanographic  Office 
(NAVOCEANO)  Digital  Bathymetric  Data  Base  -  Variable  Resolution  (DBDB-V)  product. 
These  archived  bathymetries  can  be  used  as  a  baseline.  At  most  sites  of  interest,  at  least 
one  high-resolution  survey  should  have  been  performed  to  provide  more  detailed 
bathymetry.  At  many  locations,  repeated  surveys  allow  an  assessment  of  bathymetric 
variability  which  is  crucial  to  predicting  the  likelihood  of  re-exposure  of  previously  buried 
munitions.  Models  of  dynamic  coastal  morphology  [Hanson  etal,  2003,  Jenkins  and  Inman, 
2006]  can  be  employed  to  estimate  the  variability  if  historic  measurements  are  not 
available.  For  the  preliminary  BN,  the  variability  is  specified  as  a  standard  deviation  of 
vertical  variation  (hSigma)  about  the  local  mean  depth  at  each  cross-shore  location  (see 
Section  3.5).  This  statistic  could  be  derived  from  repeated  surveys,  or  computed  from 
Monte-Carlo  exploration  of  models  of  cross-shore  profile  evolution  in  response  to  regional 
wave  climatology. 
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Determining  the  appropriate  width  for  the  state  intervals  of  the  WaterDepth  node 
requires  the  consideration  of  both  the  measurement  accuracy  of  h  and  the  sensitivity  of  the 
child  node  Ubot.  The  average  vertical  resolution  of  bathymetry  across  the  region  of  interest 
varies  according  to  its  source.  Modern  bathymetry  derived  from  airborne  lidar  should 
achieve  accuracy  of  0.15  m,  however,  substantially  more  areas  are  covered  by 
hydrographic  surveys,  where  recent  standards  for  sonar  soundings  only  require  a  vertical 
uncertainty  of  less  than  0.5  m  in  shallow  (<  100  m)  areas  [Wells  and  Monahan,  2002]. 
There  is  a  corresponding  increase  in  horizontal  resolution  available  from  coastal  lidar 
studies  (pixel  size  of  2  to  3  m2)  higher  than  that  available  from  archived  hydrographic 
sources  (on  the  order  of  100  m2).  Figure  3.4  shows  an  example  QGIS  raster  layer  of  water 
depth  formed  from  a  combination  of  hydrographic  and  lidar  sources.  The  bathymetry 
shown  illustrates  this  discrepancy  in  resolution:  the  regional  bathymetry  is  taken  from 
DBDB-V  which  has  3  second  resolution,  i.e.  pixels  every  ~93  m.  A  subsection  of  the 
nearshore  region  (the  beach  out  to  less  than  8  m  depth),  is  formed  from  lidar  data  taken  in 
2005,  in  a  USACE  Joint  Airborne  Lidar  Bathymetry  Technical  Center  of  Expertise  (JALBTCX) 
survey  after  Hurricane  Katrina,  which  has  resolution  of  close  to  2  m2,  and  vertical  accuracy 
between  15  to  30  cm  [Guenther  etal,  2000]. 


The  bottom  orbital  velocity  is  proportional  to  the  reciprocal  of  the  square  root  of  h  (given 
the  same  wave  height),  therefore  Ubot  varies  slowly  in  relatively  deep  water.  For  h  >  1  m,  it 
requires  a  decrease  of  over  2  m  to  cause  an  increase  in  Ubot  larger  than  the  state  intervals 
specified  in  Table  3.1.  In  very  deep  water,  the  width  of  the  interval  can  be  even  larger. 
However,  in  shallow  water,  the  orbital  velocity  becomes  substantially  more  sensitive  to 
changes  in  h,  so  that  the  water  depth  interval  must  be  as  small  as  0.5  m  to  resolve  Ubot 
sufficiently.  The  choice  of  state  intervals  for  the  WaterDepth  node  are  listed  in  Table  3.1 
and  illustrated  in  Figure  3.2.  Note,  again,  that  the  value  of  h  is  considered  to  be  fixed  in 
time.  Also,  for  any  given  spatial  cell,  the  value  lies  completely  within  a  single  state  interval; 
in  fact,  that  is  how  the  boundaries  of  that  spatial  cell  were  defined  (see  Section  2.2).  The 
color  scale  in  Figure  3.4  uses  the  water  depth  state  intervals,  so  that  the  color  changes 
delimit  the  isobars  where  there  is  a  state  change. 
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Figure  3.4  Water  Depth  near  Panama  City,  Florida.  Large  scale  bathymetry  from  DBDB-V. 
High-resolution  near  shoreline  from  2005  survey  by  USACE  JALBTCX.  Red  dots  mark  the 
locations  of  the  TREX13  field  experiment  locations  [Calantoni,  2015]. 


3.3  Input  Nodes:  Wave  Forcing 

The  nodes  describing  the  characteristics  of  wave  forcing  must  capture  the  temporal 
variation  of  the  dynamic  forcing  which  causes  UXO  burial  and  migration.  The  input  PDF  of 
wave  heights  should  be  computed  as  a  histogram  over  a  time  period  that  is  long  enough  to 
be  representative  of  the  conditions  at  the  site.  For  the  probabilistic  approach  taken  by 
UnMES,  the  exact  burial-motion  sequence  phasing  (Section  2.3)  is  assumed  to  be 
unimportant;  therefore,  for  adequate  statistics,  the  PDF  should  be  developed  over  a  time 
period  encompassing  multiple  cycles  of  burial-motion  as  well  as  multiple  cycles  of  wave 
patterns  (diurnal,  seasonal,  annual).  Figure  3.5  shows  an  example  histogram  of  wave  height 
computed  from  8  years  of  wave  buoy  measurements  offshore  of  FRF  at  Duck,  NC  [FRF, 
2011].  Note  at  any  moment  in  time  an  entire  spectrum  of  waves  is  present;  the  single 
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statistic  represented  here  is  significant  wave  height,  HSig,  computed  as  4 a  of  the  sea  surface 
elevations  measured  during  the  measurement  window  (usually  V2  hour). 
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Figure  3.5  Probability  distribution  of  significant  wave  height  computed  from  measurements 
at  the  FRF  Waverider  buoy:  a)  shows  the  histogram  of  wave  heights  measured  at  a  location 
in  a  water  depth  of  17  m  during  the  years  1998  to  2006  with  magenta  line  the  best  fit 
lognormal  distribution,  b)  transformed  wave  heights  at  7  m  depths;  green  line  is  the  best  fit. 
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3.3.1  Shallow  Water  Wave  Transformation 

For  most  coastal  locations,  instruments  measuring  wave  conditions  are  located  offshore  in 
deep  water.  Many  U.S  coastal  sites  are  serviced  by  offshore  wave  buoys  maintained  by 
National  Data  Buoy  Center  (NDBC),  part  of  the  National  Oceanic  and  Atmospheric 
Administration  (NOAA).  Offshore  wave  model  predictions  are  available  at  locations  in 
between  buoys  from  the  operational  wave  model  WAVEWATCH  III  (Tolman,  2009),  with 
hindcast  re-analyses  archived  back  to  1999.  To  populate  the  cells  of  the  GIS  framework  for 
UnMES,  which  is  generally  most  concerned  with  shallower  regions  close  to  shore,  the 
transformation  of  the  shoaling  waves  must  be  computed.  NOAA’s  National  Centers  for 
Environmental  Prediction  (NCEP)  is  currently  developing  the  Nearshore  Wave  Prediction 
System  (NWPS)  to  provide  high-resolution  nearshore  wave  products  [van  der  Westhuysen, 
2012].  Until  NWPS  is  routinely  operational,  the  shoaling  of  deep-water  waves  to  the 
shallow  region  of  most  interest  to  UnMES  will  require  using  a  gridded  nearshore  wave 
transformation  model  such  as  SWAN,  which  was  developed  at  Delft  University  and  is  being 
incorporated  into  NWPS  [Settelmaier  et  ah,  2011],  or  the  US  Army  Corps  of  Engineers 
(USACE)  model  STWAVE  [Smith  etal,  1999]. 

Figure  3.6  illustrates  the  effect  of  shoaling  on  two  waves  approaching  the  beach  at  FRF, 
apply  the  wave-evolution  model  by  Thornton  and  Guza  [1983].  The  water  depth  (in  brown) 
shows  a  pronounced  longshore  bar  at  200  m  offshore,  then  increases  with  a  slope  of  about 
0.008  out  to  the  depth  of  8  m  where  a  wave  gauge  is  located.  An  offshore  wave  height  of  2.4 
m  is  plotted  in  dark  blue;  this  wave  can  be  seen  to  steadily  decrease  in  height  as  it  moves 
inshore,  its  energy  dissipating  due  to  the  combined  effect  of  bed  friction  and  depth-induced 
breaking.  At  the  crest  of  the  longshore  bar,  this  wave  breaks,  losing  height  rapidly,  then 
continuing  with  HSig  about  0.75  m  to  just  before  the  shore.  In  comparison,  a  wave  in  8  m 
with  Hsig  =  0.72  m,  shown  in  cyan,  does  not  dissipate,  but  experiences  slight  HSig  growth  as 
the  shallower  depths  slow  the  wave,  reducing  its  wave  length,  causing  a  compensating 
increase  in  height.  Breaking  again  occurs  at  the  crest  of  the  longshore  bar.  Overlaid  on 
Figure  3.6  are  dashed  vertical  lines  marking  the  discretization  intervals  for  the  water  depth 
parameter  determined  in  the  Section  3.2.2.  It  can  be  seen  that  in  order  to  spatially  resolve 
the  nearshore,  closely-spaced  cells  are  required.  However,  the  near-surf-zone  is  a  region  of 
large  uncertainty,  so  that  such  high  resolution  may  not  be  warranted. 
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Figure  3.6  Cross-shore  bathymetry  and  shoaling  wave  heights  at  FRF,  Duck  NC.  Two  example 
waves  are  shown  with  offshore  heights  of  2.4  m  and  0.7  m  (marked  by  magenta  dots). 


The  use  of  Bayesian  networks  to  model  wave  evolution  from  offshore  into  the  surf  zone  has 
been  explored  [Plant  and  Holland,  2011]  and,  in  the  future,  UnMES  could  be  expanded  to 
handle  that  portion  of  the  required  modeling  effort.  For  the  current  version,  it  is  assumed 
that  wave  characteristics  at  all  spatial  cells  are  determined  offline,  and  provided  to  UnMES 
as  an  input  PDF. 

When  modeling  a  single  location,  the  prior  distribution  of  the  HSig  input  node  can  be 
approximated  using  a  histogram  such  as  shown  in  Figure  3.5.  To  provide  wave  information 
for  all  the  cells  in  a  spatial  implementation  of  UnMES,  it  is  more  practical  to  represent  the 
PDF  of  Hsig  parametrically,  by  fitting  the  observations  to  a  parametric  continuous 
probability  distribution  using  maximum  likelihood  estimation.  Then  the  GIS  input  consists 
of  one  or  two  layers  of  the  appropriate  parameter  values.  Netica  provides  an  extensive 
library  of  built-in  functions,  including  many  analytic  probability  distributions.  Wave 
formation  via  energy  cascade  is  a  multiplicative  process;  wave  heights  are  all  positive 
values  and  their  distribution  is  highly  skewed,  with  a  low  mean  value  but  a  large  variance. 
Therefore  the  lognormal  distribution  is  a  reasonable  approximation,  requiring  two 
parameters,  p,  and  o.  Overlaid  on  Figure  3.5a  in  magenta  is  the  best  fit  lognormal 
distribution  for  the  wave  heights  observed  at  the  offshore  FRF  buoy  in  17  m  depth.  The 
green  line  in  Figure  3.5b  shows  the  wave  height  distribution  after  application  of  the 
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Thornton  and  Guza  [1983]  shoaling  transformation  for  a  nearshore  location  with  a  water 
depth  of  7  m.  The  lognormal  parameterization  is  satisfactory  for  representing  the  bulk  of 
the  wave  distribution  at  this  site.  In  order  to  precisely  model  the  tail  of  the  distribution  (the 
largest  wave  heights),  which  represents  the  most  important  forces  driving  burial,  the  use  of 
an  exceedance  probability  formulation  will  be  investigated. 


3.3.2  Discretization  of  Wave  Height  and  Period 

For  fixed  water  depth  h,  bottom  orbital  velocity  is  approximately  linearly  proportional  to 
wave  height,  Ubot  -  m  *  HSig  ,  with  a  secondary  dependence  on  wave  period  T.  The  slope  m 
approaches  1  in  shallow  water  (h  <  2  m).  This  strong  sensitivity  of  the  child  node  Ubot  to 
changes  in  HSig  constrains  the  state  intervals  of  the  WaveHeight  node  to  be  fairly  small.  At 
a  given  spatial  location  in  UnMES,  h  is  considered  fixed,  but  wave  height  and  period  will 
vary  over  time.  Figure  3.7  illustrates  the  sensitivity  of  Ubot  to  Hsig  and  T  for  two  water 
depths.  In  Figure  3.7a  the  depth  is  shallow  ( h  =  2.5  m)  and  Ubot  is  seen  to  increase  rapidly 
with  HSig,  with  m  =  0.6  at  T=  3  s,  up  to  m  >  0.9  for  T  >  6  s.  To  resolve  this  relationship 
requires  very  small  state  widths  for  the  WaveHeight  node,  with  Hsig  intervals  as  small  as 
10  cm.  However,  the  maximum  wave  height  in  the  CPT  is  only  2.6  m;  all  larger  waves  are 
breaking  (plotted  as  white),  and  can  be  assigned  to  the  "surf  zone”  category.  There  is  very 
little  dependence  on  wave  period,  so  that  node  can  adequately  be  represented  by  only  two 
states,  below  and  above  T  =  4  s.  This  Hsig  discretization  for  a  shallow  water  depth  is 
illustrated  in  Figure  3.2  and  marked  by  the  black  arrows  along  the  x-axis  of  Figure  3.7a. 

In  deeper  water  with  h-  7  m  (Figure  3.7b),  Ubot  increases  more  slowly  with  Hsig;  m  ranging 
from  0.15  to  0.5.  The  maximum  Hsig  to  delimit  is  7.6  m;  all  larger  waves  would  break  before 
reaching  this  depth.  The  state  interval  boundaries  for  Hsig  are  computed  at  the  Ubot  interval 
contours  (plotted  in  white)  for  the  wave  period  (T  =  6  s),  which  is  the  mean  T  for  the 
example  Gulf  of  Mexico  location  shown  in  Figures  3.3  and  3.4.  The  variation  over  wave 
period  at  this  depth  will  be  captured  with  three  states.  Note  that  wave  period  is  invariant;  it 
does  not  change  as  the  waves  shoal  towards  shore,  and  can  therefore  be  represented  by  the 
same  distribution  at  every  cell.  It  is  clear  that,  given  local  water  depth,  the  BN  at  each 
spatial  cell  can  be  designed  to  adequately  resolve  the  relevant  wave  height  forcing  using 
the  number  of  intervals  that  were  chosen  for  the  bottom  velocity  node.  The  resulting 
choices  of  state  intervals  are  listed  in  Table  3.1  for  the  two  example  water  depths.  In  this 
manner,  the  CPTs  for  the  nodes’  relationships  can  be  represented  using  a  manageable 
number  of  entries. 
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3.3.3  Wave  Direction 

Due  to  refraction,  shoaling  waves  are  constrained  to  increasingly  be  directed  toward  the 
shoreline  (angle  z  =  0°).  Because  it  has  been  shown  in  some  locations  that  shore  erosion  is 
sensitive  to  directional  changes  in  the  incoming  waves  [Ruggiero  et  ah,  2006]  an  input  node 
for  the  incident  angle  is  included  in  UnMES.  Also,  the  longshore  current  direction  is  driven 
by  this  angle  and  used  for  predicting  the  UXO  migration  direction.  Measurements  of  wave 
direction  are  not  available  as  readily  as  those  for  height  and  period.  The  percentage  of 
waves  that  approach  obliquely  from  either  the  east  or  the  west  (assuming  a  shoreline  along 
constant  latitude)  may  be  estimated  from  beach  observations  including  video  monitoring 
[Holman  and  Stanley,  2007].  In  this  version  of  the  UnMES  the  wave  direction  parameter  is 
used  as  a  general  indicator,  therefore  a  coarse  discretization  of  3  states  is  chosen,  with  the 
state  boundaries  at  z  =  ±20°. 


3.4  Input  Nodes:  Distribution  of  UXO 
3.4.1  UXO  Type  &  Abundance 

The  exact  composition  of  munitions  at  a  site  is  usually  not  well  known.  However  an 
estimate  of  the  type  of  munitions  present  is  critical  to  the  prediction  of  the  current  state  of 
UXO  burial  and  accumulation  as  their  size,  shape,  and  density  affect  burial  and  migration 
behavior.  The  size  of  munitions  of  interest  ranges  from  small  20  mm  bullets,  through 
projectile  shells  and  mortars  with  diameters  from  40  to  81  mm,  and  howitzers  as  large  as 
155  mm.  Diameter,  length  and  weight  data  for  common  Naval  munitions  are  readily 
available,  and  several  studies  have  relied  on  estimating  density  by  assuming  a  cylindrical 
shape.  However  many  UXO  of  concern  have  more  complex  shapes.  (Access  to  additional 
munitions  details  in  DoD  restricted  ordnance  databases  will  be  beneficial  for  the 
development  of  UnMES.)  Preliminary  site  surveys  can  provide  a  prior  probability  for  the 
input  nodes  for  UXO  Type  and  UXO  Abundance.  At  some  sites,  such  as  Camp  Perry,  OH, 
there  are  historical  firing  documents  that  allow  estimates  of  the  likely  artillery  areal 
concentration  [Wilson  et  al.,  2008b].  In  this  version  of  UnMES,  UXO  Types  include  the 
category  "bullet”  with  a  diameter  of  22  mm  and  a  high  density  (specific  gravity,  Sg  =  7). 
The  type  "shell"  is  assigned  an  intermediate  diameter  of  10  cm  and  is  considered  to  have  a 
tapered  shape  and  moderate  density  (Sg  =  4.5),  while  "bomb”  represents  a  larger, 
cylindrical  UXO  with  Sg  =  2.7.  Figure  3.2  illustrates  the  UXO  Type  node  instantiated  with 
100%  of  the  type  "shell". 
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Figure  3.7  Bottom  orbital  velocity  Ubot  as  a  function  of  wave  height  Hsig  and  wave  period  at 
two  example  water  depths:  a)  h  =  2.5  m  and  b)  h  =  7  m.  The  black  arrows  mark  the  state 
interval  boundaries  required  to  discretize  Hsig  with  adequately  resolution  in  Ubot. 
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3.4.2  Initial  Burial 

The  possibility  of  burial  upon  impact  in  sandy  sediments  is  very  low  for  munitions  that 
were  dumped  at  sea  off  a  ship,  based  on  research  performed  for  the  ONR  Mine  Burial 
program  [Rennie  et  al.,  2007].  The  hardness  of  a  sand  seabed,  usually  measured  as 
sediment  shear  strength,  is  quite  high  (>12  kPa),  and  the  impact  velocity  of  ship-dropped 
munitions  is  usually  low,  less  than  10  m/s.  However,  air-dropped  bombs  or  fired  artillery 
shells  could  have  impact  velocities  exceeding  100  m/s.  Future  work  will  investigate  the 
use  of  a  proposed  model  for  high-velocity  sediment  penetration  [Chu  et  al.,  2011].  In  the 
current  version  of  UnMES,  the  Impact  Burial  node  is  a  placeholder,  indicating  minimal 
burial  (i.e.,  the  state  "less  than  10%"  instantiated). 

3.5  Input  Nodes:  Probability  of  Erosion  or  Accretion  of  Sediments 

The  node  Erosion/Accretion  in  the  BN  represents  our  knowledge  of  changes  in  the  seabed 
level  on  spatial  scales  larger  than  the  immediate  area  around  the  UXO.  There  are  several 
mechanisms  under  investigation  that  could  cause  buried  munitions  to  re-surface  or  be 
excavated.  Most  of  these  mechanisms  also  have  the  opposite  potential  and  can  bury  the 
munitions  by  covering  them  with  sediment.  At  very  long  time  scales,  the  shoreline  erosion 
is  an  obvious  mechanism  that  can  expose  previously  buried  objects.  Within  an  annual  time 
frame,  the  dominant  process  on  many  coastal  regions  is  the  shoreface  equilibrium  profile 
adjustment  to  seasonal  variation  in  wave  climate.  At  many  wave-dominated  coastal 
nearshore  sites,  the  evolution  of  surf-zone  sand  bar  position  is  the  largest  cause  of  change 
in  sediment  depth  [Plant  et  al.,  1999].  At  a  range  of  scales  from  weeks  to  years,  the 
migration  of  sand  ridges  (also  referred  to  as  sand  waves  or  megaripples)  contributes  to 
bathymetric  variation. 

It  is  a  challenge  to  deterministically  model  these  erosion/accretion  patterns  which  are 
driven  by  far-field  mechanisms  occurring  at  very  much  longer  time  scales  (months  to  years, 
or  even  decades)  than  the  elementary  processes  from  which  they  are  built,  which  occur  on 
time  scales  of  wave  events  (hours  to  days)  [Hanson  etal.,  2003].  For  the  purpose  of  UnMES, 
it  is  required  to  estimate  the  probability  that  the  amplitude  of  bathymetric  variation  will  be 
comparable  to  the  dimensions  of  the  munitions  of  concern  within  some  fairly  extended 
time  period.  To  model  long-term  migration  probabilities,  the  exact  temporal  phase  of  short 
time-scale  processes  may  not  be  required.  For  the  preliminary  UnMES,  an  empirical 
approach  is  taken,  with  erosion  or  accretion  (E/A)  specified  statistically,  assuming  a 
Gaussian  distribution  of  seabed  levels.  The  input  value  is  given  as  a  standard  deviation  of 
vertical  variation  about  the  local  mean  depth  (hSigma).  The  choice  of  the  time  period  over 
which  the  variation  is  determined  sets  the  temporal  scale  of  the  model.  Similar  to  the 
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treatment  of  dsed,  for  the  GIS  implementation  a  single  value  (hSigma  =  constant)  is 
appropriate  within  each  spatial  cell. 


Figure  3.8  shows  the  annual  variation  of  cross-shore  bathymetry  at  FRF,  derived  from  15 
years  of  survey  data  [FRF,  2011].  The  variation  about  the  mean  within  each  year,  is 
computed,  and  then  averaged  over  the  15  samples.  The  standard  deviation  hSigma  is  plotted 
in  red  over  the  mean  bathymetry  in  brown.  It  should  be  noted  that  the  standard  deviation 
computed  over  the  entire  15  year  period  is  about  twice  as  large  as  the  annual  hSigma  shown 
here,  so  the  choice  of  time  scale  is  an  important  factor  in  setting  the  probability  distribution 
of  erosion  and  accretion. 


At  two  locations,  marked  by  black  dots,  the  histograms  of  depths  are  compared  with  the 
Gaussian  PDFs  used  to  represent  them  in  the  E/A  node  (Figure  3.8b  and  c).  The  fits  could 
be  improved  by  including  a  second  parameter  quantifying  skewness,  if  the  quantity  and 
quality  of  the  bathymetry  data  warranted  the  added  complexity.  The  location  of  the 
sandbar  is  characterized  by  a  maximum  in  hSigma,  where  changes  in  seabed  level  of  over 
0.35  m  occur  relatively  frequently  (about  one  third  of  the  time).  A  change  of  this  magnitude 
is  enough  to  cover  or  uncover  many  UXO.  Farther  offshore,  hSigma  should  approach  zero, 
marking  the  depth  of  closure,  where  the  seabed  is  static  over  the  time  scale  considered.  The 
standard  deviation  computed  here  does  not  decrease  to  zero,  but  asymptotes  to  about  0.1 
m;  most  likely  that  is  the  measurement  error  for  these  bathymetric  surveys.  In  Figure  3.8a, 
the  boundaries  of  spatial  cells,  based  in  the  water  depth  discretization,  are  again  indicated 
by  the  dashed  black  lines. 

The  E/A  node  is  then  populated  with  a  normal  distribution  with  mean  zero  using  Netica’s 
internal  equation  capability.  The  E/A  node  must  be  discretized  at  intervals  fine  enough  to 
discriminate  changes  in  percentage  burial  for  the  UXO  of  interest:  5  cm  bins  are  used  out  to 
±  45  cm,  then  larger  intervals  (0.5m)  to  encompass  the  tails. 

A  separate  node,  labeled  ShorelineRetreat  (see  Figure  2.1),  is  included  in  the  BN  to 
capture  the  process  of  shoreline  erosion  and  its  contribution  to  UXO  unburial,  which  may 
be  important  at  some  sites.  For  example,  at  Camp  Perry  on  Lake  Erie,  the  shoreline  retreat 
has  been  documented  at  rates  of  about  2  m  per  year,  and  appears  to  be  a  major  contributor 
to  the  problem  of  UXO  reappearance  on  the  beach  [Wilson  etal,  2008b],  While  this  process 
could  be  represented  using  a  negative  skewness  in  the  E/A  distribution,  it  is  clearer  to 
break  shoreline  retreat  out  as  a  separate  variable.  While  the  E/A  distribution  can  be 
developed  from  data  or  model  results  covering  a  yearly  to  decadal  time  scale,  the  rate  of 
shoreline  erosion  must  usually  be  estimated  by  extrapolation  of  large-scale  data  sets 
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[Thieler  and  Hammar-Klose,  1999],  or  the  use  of  long-term  profile  evolution  models 
developed  to  simulate  the  response  to  sea  level  rise  [Hanson  et  al.,  2003].  At  some  sites, 
rapid  erosion  can  be  triggered  by  the  construction  of  structures  such  piers  or  groins. 
Estimates  from  sediment  transport  models  customized  to  the  local  beach  engineering 
should  provide  the  mean  value  to  use  at  each  spatial  cell  in  the  UnMES  GIS. 


Figure  3.8.  Example  annual  bathymetric  variability  statistics  from  a  15  year  measurement 
record  at  FRF,  Duck,  NC. 


3.6  Intermediate  Nodes:  Scour  Burial  and  Total  Burial 

The  scour  burial  node  has  two  parents,  Ubot  and  UXO  Type  (Figure  3.2).  The  state  intervals 
for  the  scour  burial  node  were  chosen  early  in  the  design  process  (Table  3.1)  and  have 
units  of  percentage,  denoting  the  percent  of  the  UXO  diameter  that  is  buried.  The  CPT  for 


28 


Johns  Hopkins 

APPLIED  PHYSICS  LABORATORY 


FPS-T- 15-0333 
Technical  Report  Page  29 


the  scour  burial  node  is  built  from  the  equations  defining  the  burial  behavior  as 
determined  by  laboratory  and  field  experiments  [Rennie  and  Brandt,  2014]  which  is 
reviewed  in  Section  4.1  and  illustrated  in  Figure  3.1.  Netica  allows  the  specification  of  a 
multipart  equation;  in  this  initial  UnMES,  the  equation  is  tripartite  with  a  different 
parameter  used  for  each  of  the  three  example  UXO  types  (Section  3.4.1).  When  Netica 
converts  an  equation  to  a  CPT,  a  decision  must  be  made  as  to  where  within  each  discretized 
state  to  evaluate  the  equation.  Netica’s  algorithm  chooses  a  number  of  points  at  random 
positions  within  the  state,  and  uses  the  average  of  the  results  obtained.  To  reduce 
sampling  uncertainty,  the  CPTs  in  UnMES  are  built  using  a  large  number,  0(1000),  of 
random  points. 


Total  burial  represents  the  interplay  between  the  predicted  scour  and  impact  burial  and 
the  erosion  or  accretion  at  the  cell.  The  parents  for  the  TotalBurial  node  are 
ScourBurial,  ImpactBurial,  E/A,  and  UXO  Type,  which  specifies  the  UXO  diameter. 
TotalBurial  is  presented  in  dimensional  units  (cm)  and  is  discretized  to  have  finer 
discrimination  at  low  burial  (state  interval  widths  of  2,  3,  or  5  cm  for  burial  less  than  20 
cm),  but  coarser  resolution  at  deeper  burial  (widths  of  10  to  25  cm  up  to  1  m  burial).  In  this 
UnMES  version,  the  equation  is  simply  the  sum  of  the  dimensional  scour  burial  and  E/A. 
Future  work  will  investigate  an  additional  parameter,  a  multiplicative  factor  applied  to 
scour  burial  to  adjust  for  the  time  scales  of  the  different  processes.  Figure  3.9  presents  an 
illustration  of  the  section  of  the  BN  with  TotalBurial  combining  the  ScourBurial  and 
E/A  distributions,  at  a  location  with  the  E/A  parameter  set  to  hSigma  =  0.12m,  and  zero 
shoreline  retreat.  The  water  depth  for  the  spatial  location  in  this  example  is  h  -  7  m,  and 
the  UXO  Type  is  "bomb”,  which  as  a  diameter  of  15  cm. 


3.7  Results  Nodes:  Migration  probability 

A  number  of  studies  have  been  conducted  to  determine  the  bottom  flow  required  to  initiate 
motion  of  sediment  particles  on  the  seabed.  The  more  relevant  literature,  involving  larger 
particles  (gravel  or  pebbles)  on  heterogeneous  sediment  mixtures,  has  been  extended  with 
laboratory  experiments  to  study  the  mobility  of  UXO-sized  objects,  as  described  in  Section 
4.2. 
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Figure  3.9  Section  of  the  UnMES  BN  showing  Burial  and  Migration  nodes.  The  wave  forcing 
(not  shown)  represents  a  strong  event  in  the  location  water  depth  of  7  m. 


3.7.1  Migration  Distance 

In  this  preliminary  UnMES  design,  the  MigrationDistance  node  is  partitioned  into 
three  states;  these  are  labeled  "Stay"  for  no  significant  movement;  "Near"  for  distances  less 
than  50  m;  and  "Far"  for  further  movement.  An  intermediate  node  called  Exposure  is 
introduced  to  determine  whether  the  UXO  is  "locked  down”  by  burial.  The  exposure  node 
has  two  states:  "exposed”  when  TotalBurial  is  less  than  half  the  UXO  diameter,  and 
"buried"  for  all  other  conditions.  The  CPT  entries  in  the  MigrationDistance  node  with 
the  parent  state  "buried"  are  all  set  to  "Stay".  If  "exposed",  then  the  probability  of  motion  is 
based  on  the  bottom  velocity  distribution  and  UXO  density  (formulated  as  the  Object 
Shields  parameter  discussed  in  Section  4.2)  compared  to  the  UXO  size  versus  the  bottom 
roughness  scale,  represented  by  the  node  DoverK. 

While  the  equation  to  compute  a  threshold  of  mobilization  is  given  in  Section  4.2,  there  is 
little  guidance  yet  on  how  far  the  UXO  will  travel  once  motion  is  initiated.  For  the  current 
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UnMES  implementation,  an  ad  hoc  formulation  is  utilized.  When  the  mobility  threshold  is 
just  slightly  exceeded  (i.e.  Object  Shields  parameter  within  an  order  of  magnitude  of  the 
threshold)  UXO  are  envisioned  to  slide  or  roll  along  the  seabed.  Laboratory  studies  of 
rolling  cylinders  [Davis  et  al.,  2007]  show  the  object  velocity  to  be  about  70%  of  the  bottom 
flow  speed.  Therefore,  for  these  cases  a  rule-of-thumb  distance  is  estimated  based  on  the 
speed  0.7-Ubot,  acting  over  half  the  wave  period.  When  the  forcing  strongly  exceeds  the 
mobility  threshold  [by  more  than  one  order  of  magnitude),  the  object  travels  suspended  in 
the  current,  and  a  full  100%  of  Ubot  is  used.  When  the  forcing  dramatically  exceeds  the 
threshold  (more  than  two  orders  of  magnitude),  it  is  assumed  that  conditions  are  such  that 
the  entire  seabed  is  mobilized,  with  burial  unlikely,  and  multiple  instances  of  mobilization 
occur.  Only  these  extreme  conditions  can  generate  migration  in  the  "Far"  state.  These  ad 
hoc  rules  make  the  MigrationDistance  node  conditionally  dependent  on  at  least  five 
parent  nodes.  To  refine  these  coarsely-determined  migration  estimates  will  require 
further  research,  especially  the  inclusion  of  results  from  recent  field  experiments  in 
energetic  conditions. 


3.7.2  Migration  Direction 

Because  generally  detailed  knowledge  is  lacking  about  the  direction  that  a  UXO  will  move, 
the  MigrationDirection  node  is  discretized  into  a  small  number  of  states. 
MigrationDirection  is  set  up  representing  directional  quadrants:  onshore  (±45°), 
"east"  (to  the  left  facing  onshore),  "west",  and  offshore  (which  is  separated  into  two  states 
to  encompass  -180  to  -135°  and  +135  to  +180°).  A  straightforward  prediction  of  migration 
direction,  proceeding  from  a  simple  assumption  of  mobilization  during  the  passage  of  a 
wave’s  crest,  will  be  onshore,  with  some  component  of  alongshore,  based  on  the  input 
wave  direction.  However,  consideration  of  the  general  tendency  of  objects  to  migration 
downslope  increases  the  likelihood  of  offshore  movement.  A  detailed  CFD  modeling  study 
using  the  UXO-MM  code  [Jenkins  et  al.,  2013]  at  a  site  Adak,  AK,  driven  by  a  20-year  wave 
record,  showed  that  the  dominant  direction  of  movement  was  offshore,  with  less  than  10% 
of  the  modeled  migration  towards  the  beach,  and  about  a  third  traveling  alongshore.  Field 
studies  over  3  to  5  months  offshore  of  FRF  at  Duck,  NC  [Wilson  et  al.,  2008a]  reported 
onshore  movement  for  12%  to  20%  of  the  UXO  surrogates,  with  offshore  migration  in  a 
third  to  a  half  of  the  observations.  Our  ability  to  predict  the  dominant  direction  of  UXO 
migration  will  require  additional  site-specific  environmental  knowledge,  such  as 
prevalence  of  the  longshore  current.  For  the  current  UnMES  implementation,  a  prior  PMF 
is  assigned  to  the  MigrationDirection  node  that  reflects  the  long-term  Adak  model 
results,  and  is  shown  in  Figure  3.9. 
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4  Core  Physics:  process  models  in  UnMES 

Process  models  for  UXO  scour  burial,  mobility  and  re-exposure  are  at  the  core  of  the 
UnMES.  Significant  progress  has  been  made  on  development  of  appropriate  models  for 
UXO  burial  and  mobility  based  on  the  extant  literature,  available  field  data  and  experiments 
designed  specifically  for  the  properties  of  UXO  (size  and  density).  A  brief  summary  of  these 
models  is  presented  below  while  the  details  of  their  development  are  presented  in  Rennie 
and  Brandt  [2014],  Appropriate  process  models  for  UXO  re-exposure  are  currently  under 
development. 

4.1  Scour  Burial  Model 

Water  flow,  either  the  oscillatory  motion  of  waves,  or  steady  currents,  is  accelerated  in  the 
vicinity  of  an  object  lying  on  the  seabed,  causing  sediments  to  erode  at  the  ends  of  the  UXO. 
Scour  pits  are  thus  formed,  and,  as  scouring  continues,  the  UXO  falls  into  the  enlarging  pit. 
The  pit  then  becomes  a  sediment  trap  and  can  fill  in.  This  process  and  the  following  process 
model  results  are  described  in  Rennie  and  Brandt  [2014], 

As  the  extant  data  on  scour  around  finite  objects  is  quite  limited.  Available  studies  include  a 
laboratory  study  by  Demir  and  Garcia  [2007]  on  cylinders  in  oscillatory  flow,  field  studies 
of  large  sea-mine  surrogates  [Trembanis  et  al,  2007].  and  a  series  of  experiments  utilizing 
UXO  surrogate  shapes  conducted  at  jHU/APL  [Rennie  and  Brandt,  2014].  In  these  recent 
studies  large  "shell"  like  cylinders,  tapered  cylinders  and  small  "bullet”  sized  cylinders 
were  utilized. 

The  results  and  prior  data  are  shown  in  Figure  4.1,  where  the  fraction  of  final  scour  burial 
is  plotted  as  a  function  of  the  sediment  Shields  parameter,  0,  which  is  the  dimensionless 
parameter  that  encompasses  the  driving  parameter  for  the  scour  process.  0  is  defined  as 


9  (.^sed  U^sed 


where  U  the  bottom  current,  g  gravitational  acceleration,  Ssed  =  ps/pw,  with  ps  the  sediment 
grain  density,  pw  the  water  density,  and  dsed  the  median  sand  grain  size.  The  friction  drag 
force  coefficient  fT  represents  the  skin-friction  shear  stress  acting  on  the  bed,  rather  than 
the  form  drag  acting  on  any  individual  sand  particle  (Soulsby,  1997). 
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Figure  4.1  shows  the  previously  compiled  laboratory  data  [Demir  and  Garcia,  2007], 
labeled  "DG2007"  and  field  data  [Trembanis  et  al.,  2007],  labeled  "ONR  MBP  FIELD."  Red 
circles  indicate  ]HU/APL  laboratory  results  using  cylinder  with  diameter  D  =  10  cm.  Brown 
triangles  show  the  burial  observed  with  the  tapered  cylinder  (D  =  7.94  cm).  The  small 
diameter  (D  =  2.54  cm)  cylinders  shown  as  magenta  circles.  The  considerable  scatter 
apparent  is  the  result  of  the  natural  variability  of  the  turbulent  scouring  process,  the 
variability  in  the  scour  pit  growth  and  differences  in  the  local  sediment  environment. 
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Figure  4.1.  Observations  of  scour  burial  of  UXO  and  sea-mine  surrogates.  Burial  is  reported 
as  the  fraction  of  the  surrogate  diameter  and  plotted  versus  the  Sediment  Shields  number  0. 

It  was  found  that  the  burial  fraction  can  be  reasonably  represented  by  the  relationship 


—  =  a  9b, 

D 


(4.2) 


where  Beq  is  the  equilibrium  burial  depth  and  a  and  b  are  empirically  determined 
constants.  This  relationship  gives  a  reasonably  simple  process  model  for  burial  depth  for 
use  in  the  UnMES  that  relates  the  UXO  properties  to  the  environmental  conditions  through 
the  Shields  parameter. 
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The  blue,  black  and  green  dashed  lines  show  the  empirical  fits  to  Equation  4.2  from 
Friedrichs  [2014],  Demir  and  Garcia  [2007]  and  Whitehouse  [1998],  respectively.  For  large, 
cylindrical  objects,  b  lies  between  0.8  and  1.0  and  a  is  about  2.2.  The  scale  factor  a  is  larger 
for  tapered  shapes  or  very  small  UXO  such  as  bullets,  which  exhibited  a  tendency  to  bury 
more  easily.  The  magenta  line  represents  the  best  fit  of  Equation  4.2  for  the  combined 
tapered  and  small  UXO  observations  with  values  of  a  -  5.9  and  b  =  1.  These  empirical 
parameters  are  applied  in  the  Netica  equations  defining  the  CPT  for  the  ScourBurial 
node  (Section  3). 


4.2  Mobility  Model 


Motion  of  a  bottom-sitting  object  results  when  the  force  on  the  object  due  to  its  drag  in  the 
local  flow  exceeds  the  resisting  force  due  to  the  objects  weight  and  inertia.  For  large 
objects  compared  to  the  bed  sand  grain  size,  the  appropriate  parameter  governing  the 
criteria  for  the  onset  motion  is  the  Object  Shields  number  0Obj  [Friedrichs,  2015]  defined  in 
analogy  to  the  Shields  parameter,  0,  as 


@  nh  i 


°bj  9  ( Sobj  ~  l)  D 


(4.3) 


where  the  object  diameter,  D,  replaces  the  characteristic  sand  grain  size,  dsed,  and  the 
specific  gravity  of  the  object,  S0bj  ,  replaces  the  specific  gravity  of  the  sand  grains  used  in 
Equation  4.1. 

Figure  4.2  shows  the  mobility  threshold,  0Obj,  data  plotted  versus  the  ratio  of  the  object 
diameter  D  to  bed  roughness,  k.  This  ratio  D/k,  parameterizes  geometrical  aspects  of  the 
force  balance,  and  is  represented  by  the  BN  node  DoverK  in  UnMES.  The  prior  extant  data 
considered  large  objects  in  field  studies  shown  in  green,  magenta  and  brown  on  the  right 
and  the  motion  of  the  particles  from  a  it  sand  bed  itself  in  black  on  the  left  (see  Rennie  and 
Brandt  [2014]  for  details).  As  UXO  fall  in  the  intermediate  range  that  was  not  covered,  a 
series  of  laboratory  studies  were  performed  to  provide  these  data,  shown  as  red  symbols. 

The  solid  blue  in  Figure  4.2  line  is  the  best  fit  power  law  for  all  data,  while  the  red  line 
shows  the  fit  for  D/k  >  4  excluding  field  data  of  natural  sediments,  with  dashed  lines 
indicating  the  95%  confidence  limit.  This  latter  curve  provides  the  process  model  needed 
in  UnMES  as 
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0Obj  =  1.2  (D/k)"0'62 


(4.4) 


Equation  4.4  is  used  in  the  definition  the  CPT  of  the  MigrationDistance  node  (Section 
3.7).  It  should  be  noted  that  the  onset  of  UXO  motion  for  objects  on  a  sand  bed  is  somewhat 
more  complicated  due  to  the  fact  that  some  scour  may  occur  prior  to  the  onset  of  motion 
which  will  alter  the  flow  pattern  around  the  UXO  and  thus  the  appropriate  value  of  0Obj,  as 
discussed  in  Rennie  and  Brandt  [2014]. 


Figure  4.2.  Compiled  data  for  the  mobility  threshold  of  objects  on  the  seabed.  JHU/APL 
laboratory  results  shown  are  in  red.  The  solid  blue  line  is  the  best  fit  power  law  for  all 
previous  and  new  data  combined;  the  red  line  shows  the  fit  for  D/k  >  4,  which  excludes  field 
data  of  near-homogenous  sediments.  The  dashed  lines  indicate  the  95%  confidence  limit 
around  the  large  D/k  fit.  The  grey  lines  indicate  the  state  interval  boundaries  for  the  Object 
Shields  parameter,  delimiting  different  defined  modes  of  migration  behavior  (Section  3.7.1). 
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4.3  Re-exposure  Models 

Predicting  the  probability  of  re-exposure  is  a  particularly  important  factor  in  estimating 
migration  potential,  as  UXO  can  be  mobilized  only  when  buried  less  than  about  half  of  its 
casing  diameter.  In  addition,  estimation  of  the  percentage  of  time  spent  covered  by 
sediment  is  important  to  predicting  the  condition  of  the  munitions’  casings.  Complete 
burial  in  anoxic  conditions  can  result  in  casings  displaying  very  little  corrosion  even  after 
decades,  whereas  prolonged  periods  of  exposure  on  the  seabed  results  in  corroded,  leaking 
UXO.  Prediction  of  the  casings’  condition  could  be  included  in  future  versions  of  UnMES 
that  tackle  the  problem  of  estimating  seawater  pollution  from  dissolved  chemical 
components  of  munitions. 

Re-exposure  is  the  result  of  temporal  changes  in  far-field  sediment  morphology  with 
different  processes  dominant  at  different  time-scales.  At  the  MR  Program  Review 
Workshop  [SERDP,  2014],  it  was  noted  that  progress  is  needed  in  modeling  far-field 
processes  in  order  to  understand  long-term  UXO  migration  patterns.  The  proposed 
continuation  of  the  present  study  [Rennie  and  Brandt  2015a],  in  collaboration  with  the 
upcoming  SERDP  Project  MR01-026  [Friedrichs,  2015]  will  develop  improved 
parameterized  model  formulations  describing  far-field  phenomenology  including  ripple 
and  dune  migration  as  well  as  larger-scale  shoreward  movement  of  the  shoreline  and 
associated  shoreface,  both  cyclical  (e.g.,  winter/summer  beach  profile  adjustments)  or  uni¬ 
directional  (e.g.,  net  coastal  erosion).  The  current  version  of  UnMES  uses  a  data-derived 
statistical  representation  of  exposure  probability  (Section  3.5)  as  a  placeholder  until  these 
far-field  process  models  are  sufficiently  developed. 

5  Example  Implementation 

Environmental  data  is  being  gathered  in  order  to  populate  the  input  GIS  layers  for  two 
example  implementations  of  UnMES:  the  Gulf  of  Mexico  (GoM)  off  the  coast  of  Panama  City 
(shown  in  Figures  3.3  and  3.4);  and  at  the  FRF  in  Duck  NC  (Figures  3.5  and  3.6).  These  sites 
are  chosen  because  field  test  observations  will  be  available  for  future  comparison  with 
UnMES  predictions  for  verification  of  the  expert  system  utility.  The  behavior  of  a  number  of 
surrogate  munitions  deployed  in  the  GoM  was  studied  during  the  TREX13  field  test 
[Calantoni,  2015]  in  April  through  May  2013.  An  earlier  field  test  was  conducted  at  FRF 
over  a  more  extended  period  of  time  in  2005  to  2006  [Wilson  etal,  2008a]. 
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In  Figures  3.2  and  3.9  sections  of  the  UnMES  core  BN  are  shown  for  the  GoM  example. 
Figure  3.2  illustrates  the  Scour  Burial  response  for  a  wave  height  distribution 
representative  of  the  coast  off  Panama  City  in  shallow  water  ( h  is  between  2.5  and  3  m). 
The  characteristics  of  each  node  are  displayed  as  "belief  bars”,  the  probability  distribution 
approximated  by  a  PMF  (normalized  histogram).  At  this  water  depth,  which  is  inshore  of 
the  sandbar  dissipation  has  reduced  all  larger  offshore  waves  so  that  HSig  <  1.4  m,  resulting 
in  a  local  peak  in  the  wave  height  distribution  in  the  state  "1.2  to  1.4  m".  The  wave  period 
distribution  is  based  on  observations,  so  that  its  relationship  with  wave  height  is 
represented.  Here  the  mean  period  is  5.9  ±  2.7  s,  which  is  displayed  at  the  bottom  of  the 
Period  node.  The  distribution  of  bottom  orbital  velocity  is  computed  by  Netica  in  the  Ubot 
node  using  the  linear  wave  equations.  The  UXO  Type  input  node  is  set  to  all  of  state  "shell" 
(Section  3.4.1).  The  seabed  sediment  is  composed  of  medium  sand  grain  size  dsed  =  0.45 
mm,  so  that  scour  process  falls  in  between  that  of  very  fine  and  coarse  sand  behavior  as 
shown  in  Figure  3.1.  Note  that  because  dsed  is  treated  as  a  Netica  constant,  that  node’s 
value  is  available  to  all  other  nodes  in  the  BN  without  a  link  (arrow). 


Because  the  Wave  Height  node  here  represents  the  distribution  of  wave  conditions  over  an 
extended  temporal  period,  the  ScourBurial  PMF  can  be  interpreted  as  the  percentage  of 
time  for  which  conditions  are  such  that  a  proud  UXO  would  become  buried  The 
ScourBurial  node  computes  a  ~27%  chance  of  very  slight  burial  (less  than  10%  buried) 
that  corresponds  approximately  to  the  combined  probability  of  the  two  lowest  states  in  the 
Ubot  distribution,  representing  bottom  flow  speeds  less  than  20  cm/s,  and  consistent  with 
Figure  3.1.  Note  that  state  intervals  for  a  node  may  have  different  widths  -  this  can  result  in 
a  PMF’s  belief  bars  appearing  oddly  irregular.  For  example,  the  first  two  ScourBurial 
states  encompass  10  percentiles,  while  the  central  two  states  are  wider,  and  therefore 
contain  a  larger  portion  of  the  PDF,  even  though  the  probability  of  degree  of  burial  is 
actually  decreasing  monotonically.  The  state  "fully  buried"  is  populated  with  all 
possibilities  where  wave  forcing  is  strong  enough  to  arrive  at  equilibrium  burial,  which  is 
quantified  as  a  burial  depth  of  115%  the  object  diameter  [Whitehouse,  1998],  Given  the 
empirically-determined  parameters  for  the  burial  behavior  of  a  "shell"  type  UXO  (Section 
4.1),  the  "fully  buried"  state  occurs  for  all  Ubot  >  60  cm/s,  and  therefore  contains  the  sum  of 
the  probabilities  over  the  six  highest  Ubot  states. 


The  section  of  the  UnMES  BN  addressing  Total  Burial  and  Migration  is  shown  in  Figure  3.9. 
In  this  example  the  wave  conditions  represent  a  storm  event  (HSig  ~  2.5  m,  WaveHeight 
node  not  shown),  in  a  water  depth  of  ~7  m.  The  UXO  Type  modeled  is  "bomb”  which  is 
larger  and  less  dense  than  the  other  types.  Initial  impact  makes  no  contribution  to  burial  in 
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this  example.  After  storm  conditions,  the  ScourBurial  node  reports  the  large  majority 
(88%)  of  the  bomb  casings  would  likely  be  more  than  half  buried.  However,  the  possibility 
of  seabed  erosion  at  the  site  reduces  that  likelihood  somewhat,  so  that  "buried"  state  in  the 
Exposure  node  is  less  than  66%.  Of  the  exposed  UXO,  most  are  predicted  to  remain  in  their 
original  positions  ("Stay"  is  the  most  probable  state  in  MigrationDistance),  while 
about  one  third  are  predicted  to  move  a  short  distance.  There  is  an  almost  zero  probability 
predicted  for  substantial  "Far"  migration.  A  comparison  of  UnMES  predictions  with  field 
observations  is  currently  in  progress  and  will  be  reported  in  the  forthcoming 
Demonstration  Report  on  the  prototype  UnMES  implementation  [Rennie  and  Brandt, 
2015b], 


6  Summary  and  Future  Work 


The  preliminary  design  for  the  Underwater  Munitions  Expert  System  (UnMES)  is  described. 
UnMES  is  a  computer-based  probabilistic  expert  system  to  predict  the  likelihood  of  burial 
and  migration  of  abandoned  underwater  munitions.  This  approach  utilizes  simple 
deterministic  models  in  a  probabilistic  framework,  combined  with  statistical 
parameterizations  to  define  the  joint  conditional  probability  space  for  UXO  behavior. 
Probabilistic  modeling  is  appropriate  given  the  inherent  uncertainties  in  a  munitions  site’s 
history  and  imperfect  knowledge  of  the  environment  and  complex  physical  processes 
governing  UXO  burial  and  mobility.  UnMES  will  be  a  component  of  a  larger  risk  assessment 
tool  developed  for  use  by  site  managers  as  remediation  guidance.  This  larger  framework 
will  require  GIS  support  for  environmental  input  and  visualization  of  spatially  varying 
assessments.  The  Bayesian  network  for  the  prototype  UnMES  makes  use  of  the  GeoNetica 
software  product,  and  a  QGIS  framework  to  interpret  and  display  map-based  information. 


The  component  variables  of  the  Bayesian  network  are  specified  as  belonging  to  input, 
output,  and  intermediate  node  sets.  The  interrelationships  between  nodes  determines  the 
causal  structure  of  the  predictive  system.  Methods  are  developed  to  determine  the 
required  discretization  of  the  continuous  variables  delimiting  the  states  of  each  node, 
balancing  the  need  for  adequate  resolution  with  operational  tractability.  Future  work  will 
investigate  the  application  of  alternative  probabilistic  programming  environments  which 
can  handle  continuous  distributions  with  techniques  for  performing  the  Bayesian  inference 
parametrically.  This  would  eliminate  the  requirement  to  make  prior  selections  of  the  state 
intervals  for  discretizing  the  PDFs. 

Environmental  data  from  two  sites  of  interest  are  compiled  in  QGIS  layers  and  are  being 
used  as  input  to  the  prototype  UnMES  implementation.  The  task  of  demonstration  and 
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verification  of  UnMES  predictions  is  being  undertaken  using  field  test  observations  from 
recent  SERDP-funded  experiments  conducted  in  wave-dominated  coastal  locations  off 
Florida  and  North  Carolina.  Additional  environmental  data  will  be  acquired  from  the  field 
experiments  near  Martha’s  Vineyard  [Traykovski,  2015],  Access  to  details  of  these  UXO 
burial  and  migration  observations  should  be  available  in  the  near  future. 


The  Underwater  Munitions  Expert  System  will  be  part  of  a  decision  support  tool  providing 
graphical  mapping  and  visualization  of  UXO  distributions  and  burial  state  for  site 
remediation  planning.  A  key  aspect  of  the  continued  development  of  UnMES  will  be 
forming  the  metrics  by  which  prediction  uncertainty  can  be  presented  and  assessed  for  the 
evaluation  of  management  alternatives.  A  QGIS  plug-in,  the  Probabilistic  Map  Algebra  Tool 
(PMAT)  [Landuyt  et  al,  2015],  is  under  evaluation  for  possible  use.  PMAT  produces 
different  formats  of  Bayesian  belief  network  maps  to  aid  in  understanding  risk,  including 
(1)  the  most  probable  state  accompanied  by  a  map  of  the  probability  of  that  reported  state; 
and  (2)  a  map  of  expected  value  along  with  a  second  map  of  standard  deviation  cumulative 
probability  below  or  above  a  specified  threshold.  These  formats  will  be  evaluated  during 
the  prototype  UnMES  evaluation. 


Use  of  the  preliminary  UnMES  design  has  made  clear  that  additional  work  needed  on 
understanding  re-exposure  probabilities.  Further  research  is  proposed  [Rennie  and  Brandt, 
2015b,  Friedrichs,  2015],  to  develop  rational  parameterized  models  for  UXO  re-exposure 
due  to  far  field  effects  such  as  seasonal  beach  profile  shaping  as  well  as  sand  wave  or  large 
bedform  migration.  Other  topics  requiring  additional  research  are  the  potential  for  impact 
burial  and  a  better  understanding  of  the  distance  that  munitions  can  travel  once  mobilized. 
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